Weighted Visit Mean and SD

library(tidyverse)
library(censusapi)
library(tigris)
library(sf)
library(mapview)
library(usmap)
library(rjson)
library(censusapi)
library(leaflet)
library(stringdist)
library(measurements)
library(units)
library(pracma)
library(plotly)
library(ggplot2)
library(radiant.data)

Sys.setenv(CENSUS_KEY="b88495f8315f45b2d100dee9ba8f4c489a7371c2")

options(
  tigris_class = "sf",
  tigris_use_cache = TRUE
)



visits_200309 <- readRDS("P:/SFBI/Restricted Data Library/Safegraph/covid19analysis/weekly-patterns/v2/sj_visits_daily_sum_hourly_03-09-20_05-24-20.rds")

visits_200525 <- readRDS("P:/SFBI/Restricted Data Library/Safegraph/covid19analysis/weekly-patterns/v2/sj_visits_daily_sum_hourly_05-25-20_05-31-20.rds")

visit_data <-
  visits_200309 %>%
  rbind(visits_200525)

scc_population <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*",
    regionin = "state:06+county:085",
    vars = c("group(B01003)")
  ) %>%
  mutate(
    origin_census_block_group =
      paste0(state,county,tract,block_group)
  ) %>% 
  rename(
    total_pop = "B01003_001E"
  ) %>%
  dplyr::select(total_pop, origin_census_block_group)


sj_blockgroups <- readRDS("sj_blockgroups.rds")


combo_visit_pop <-
  visit_data %>%
  left_join(
    scc_population %>%
      filter(origin_census_block_group %in% sj_blockgroups$origin_census_block_group),
    by = "origin_census_block_group"
  )

date_pull <-
  combo_visit_pop %>%
  pull(date) %>%
  unique()

day_weighted <- 
    1:length(date_pull) %>% 
    map_dfr(function(i){
      
      date_filter <-
        combo_visit_pop %>%
        filter(date == date_pull[i]) 
      
      single_day_process <-
        date_filter %>%
        mutate(
          multiplied = total_visits_avg * total_pop
        ) %>%
        group_by(date) %>%
        summarise(sum_multiplied = sum(multiplied), sum_total_pop = sum(total_pop))
      
      single_day_mean <-
        single_day_process %>%
        mutate(
          weighted_visit_avg = sum_multiplied / sum_total_pop,
          mean_percapita = weighted_visit_avg / mean(date_filter$total_pop)
          
          ) %>%
        dplyr::select(date,weighted_visit_avg,mean_percapita)
      
      single_day_weights <-
        date_filter %>%
        mutate(
          weight_value = total_pop / sum(date_filter$total_pop)
        )
      
      single_day_all <-
        single_day_mean %>%
        mutate(
          weighted_visit_sd = weighted.sd(date_filter$total_visits_avg,single_day_weights$weight_value, na.rm = TRUE),
          sd_percapita = weighted_visit_sd / mean(date_filter$total_pop)
        )
        
    })

block_group_visits <-
  combo_visit_pop %>%
  mutate(
    percapita = total_visits_avg / total_pop
  ) %>%
  dplyr::select(origin_census_block_group,date,percapita)

weekends <-
  combo_visit_pop %>% 
  filter(!duplicated(date)) %>% 
  arrange(date) %>% 
  mutate(
    x = 
      case_when(
        (date %>% as.numeric()) %% 7 == 1 ~ date + 1,
        (date %>% as.numeric()) %% 7 == 4 ~ date - 1,
        TRUE ~ date
      ),
    y = 
      ifelse(
        (date %>% as.numeric()) %% 7 %in% c(2,3),
        10,
        0
      )
  ) %>% 
  dplyr::select(x,y)

weight_visit_plot <-
plot_ly(
    data = day_weighted,
    x = ~date
  ) %>%
  add_trace(
      x = weekends$x,
      y = weekends$y,
      type = "scatter",
      mode = "lines",
      fill = "tozeroy",
      fillcolor = "rgba(112,112,112,0.25)",
      line = list(color = "transparent"),
      name = "Weekends",
      showlegend=T
    )%>%
    add_trace(
      y = ~mean_percapita-sd_percapita,
      type = "scatter",
      mode = "lines",
      line = list(color = "transparent"),
      name = "City, +/- 1 S.D.",
      showlegend=F
    ) %>% 
    add_trace(
      y = ~mean_percapita+sd_percapita,
      type = "scatter",
      mode = "lines",
      fill = "tonexty",
      fillcolor = "rgba(31,164,99,0.25)",
      line = list(color = "transparent"),
      name = "City, +/- 1 S.D.",
      showlegend=T
    ) %>% 
    add_trace(
      y = ~mean_percapita,
      type = "scatter",
      mode = "lines",
      line = list(color = "rgb(0,0,0)"),
      name = "City, Weighted Average",
      showlegend=T
    ) %>%
    layout(
      margin = list(t = 25),
      paper_bgcolor='rgb(255,255,255)', 
      plot_bgcolor='rgb(229,229,229)',
      xaxis = list(
        title = "",
        gridcolor = 'rgb(255,255,255)',
        showgrid = TRUE,
        showline = FALSE,
        showticklabels = TRUE,
        tickcolor = 'rgb(127,127,127)',
        ticks = 'outside',
        zeroline = FALSE,
        fixedrange = T,
        range = c(min(combo_visit_pop$date),max(combo_visit_pop$date)+1)
      ),
      yaxis = list(
        title = "Mean Visits per Capita",
        gridcolor = 'rgb(255,255,255)',
        showgrid = TRUE,
        showline = FALSE,
        showticklabels = TRUE,
        tickcolor = 'rgb(127,127,127)',
        ticks = 'outside',
        zeroline = FALSE,
        fixedrange = T,
        range = c(0,2)
      ))

weight_visit_plot